S-wave superconductivity near a surface 
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We study the superconducting order parameter near a sur- 
face with the Bogoliubov-de Gennes formalism. For definite- 
ness we use the attractive Hubbard model. Near a surface, the 
order parameter and the density distribution exhibit "Friedel- 
like" oscillations. Although the local density of states is quite 
different from that in the bulk, the energy gap in the spec- 
trum on a surface is almost the same as the bulk value. In the 
low-density limit, however, the energy gap tends to vanish on 
a surface. 



I. INTRODUCTION 

Experimental probes of various properties of supercon- 
ductors can generally be classified into two categories - 
bulk and surface probes. Tunneling and photoelectron 
spectroscopy (PES) are examples in the second category, 
and while PES has only recently been utilized heavily to 
investigate superconductivity [1] , tunneling has long been 
an important tool for determining the energy gap [2] and 
electron-phonon related structure [3] in superconductors. 
Through comparison with results from bulk probes, it has 
been well established that the superconducting order pa- 
rameter in conventional superconductors remains robust 
at surfaces. 

The purpose of this paper is to investigate the extent 
to which this statement is theoretically understood. In 
order to do so, in the context of conventional s-wave su- 
perconductivity, we adopt a minimal model to describe 
s-wave superconductivity, namely the attractive Hubbard 
model. This choice allows us to explore various aspects 
of the problem that otherwise might be difficult to study: 
electron density dependence and weak-to-strong coupling 
tendencies, to name a few. The tight-binding and lat- 
tice features of the this model also make it amenable 
to numerical treatment. The presence of surfaces can 
be included in this model simply by requiring that elec- 
trons cannot hop beyond the surfaces, e.g., by imposing 
open boundary conditions (OBC) as opposed to periodic 
boundary conditions (PBC). 

The paper is organized as follows. In the next section 
we outline the formalism used in this study. We simply 
solve the Bogoliubov-de Gennes (BdG) equations [4], as 



described in a previous publication [5] . We also compare 
results with those obtained from the Anderson prescrip- 
tion [6,5], for which a brief review is provided. In the 
following section we present results for the order param- 
eter as a function of distance from a surface for various 
parameter regimes. It is clear that a large spatial varia- 
tion occurs, dependent on coupling strength and electron 
density. We also examine the local density of states as a 
function of distance from the surface, and find significant 
variation close to the surface. In this work we focus our 
study on one dimension, in which the impact of a surface 
can be visualized clearly. The effects of dimensionality 
on the "Friedel-like" oscillations are discussed, and also 
some results for the local density of states are shown for 
a two-dimensional system. 



II. FORMULATION 

The use of the Bogoliubov-de Gennes (BdG) equations 
[4,8] (or a corresponding quasi-classical formulation [9]) 
is required to solve for the superconducting order pa- 
rameter in cases where the order parameter varies with 
position. This will generally occur in the vicinity of a 
defect or impurity, but first and foremost, will always oc- 
cur near a surface. Even if the surface does not modify 
the material parameters, its existence causes an abrupt 
change in the effective pair potential, and hence the order 
parameter will vary in its vicinity. 

What governs the length scale of the variation in the 
order parameter as one moves away from the surface? 
Three length scales enter the problem (we leave out 
the actual thickness of the sample, which enters in the 
Tomasch-Rowell experiment [10], or in nanograin super- 
conductors). These are (i) the lattice spacing, a, (ii) 
the coherence length, £qj an d (hi) a measure of the 
inter-electron spacing, say, hp 1 , where kp is the Fermi 
wavevector. Near half-filling (i) and (iii) coincide; in the 
low-density limit, we will find that (iii) governs the or- 
der parameter variation. In all cases, one can observe 
a 'healing length', over which the order parameter at- 
tains its bulk value. This is qualitatively described by 
the coherence length, determined here with the BCS ex- 
pression: 
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where (R 2 ) is the mean square radius of an electron pair 
[11,12]. Strictly speaking, as defined this length has little 
to do with coherence (defined for a single pair). However, 
at low temperatures, to which we confine this study, this 
tends to coincide with the Ginzburg-Landau coherence 
length, which more properly describes the length scale of 
fluctuations in the order parameter away from its mean 
value. 

The BdG equations are [8,5]: 

E t ui{i)= ^A u , Ul (i') + Vm(i) + AiVi(i) (2) 

v 

Emit) = - ^AuMi') -ViVi(i) + A* Ul (i) (3) 



the non-interacting problem. For tight-binding systems 
with open boundary conditions, these solutions can be 
found analytically. Once the single-particle states are ob- 
tained, the BCS-like equations can be solved much more 
easily than their BdG counterparts. As shown previously 
[5] and further in this work, they provide an excellent 
qualitative description of the position dependence of the 
order parameter and the density of states, for a fraction 
of the numerical cost for BdG. 



III. RESULTS 



A. BdG vs. Anderson 



where 



(4) 



The self-consistent potentials, Vj and Aj, are given by 
A i = \U\Y,ui{i)v* l {i)(l-2f l ) (5) 



Vi 



l«i(i)l 2 /i + l«j(OI 2 (l-/i) 



(6) 



The index I is used to label the eigenvalues (there are 
2N of them), the index i to label the sites (1 through 



N), and the composite eigenvector is given by 



of 



total length 2N. The sums in Eqs. (5,6) are over positive 
eigenvalues only. The fi is the Fermi function, with ar- 
gument pEi, where (3 = l/ksT with T the temperature, 
and ks the Boltzmann constant. The single-site electron 
density rij, is given through Eq. (6) by Vi — —\U\rii/2. 
As usual, the parameters t and \U\ are the hopping inte- 
gral and attractive coupling strength, \i is the chemical 
potential, and the subscript in Eq. (4) refers to near- 
est neighbour sites in units of the lattice spacing. 

These equations are solved by iteration until full self- 
consistency is achieved. If (as is often the case) a given 
electron density is desired, a second iteration loop is in- 
serted to vary the chemical potential fi, until the required 
average density n = jj rij is attained. The quasipar- 
ticle amplitudes are used to calculate the local density of 
states (LDOS) given by 



Ai(cj) = J2 \ M*)| 2 6(w - Ei) + h«| 2 5(u> + Ei) 



(7) 



We will make reference to solutions to the Anderson 
equations in what follows. As described in the original 
paper [6] and explained for the attractive Hubbard model 
in Rcf . [5] , these are a BCS-like set of equations for an in- 
homogeneous system, for which solutions are required for 



In the absence of a surface or impurity, momentum 
is a good quantum number, and the BCS gap equations 
[13] can be solved. For the attractive Hubbard model, 
the resulting order parameter turns out to be constant in 
momentum space as well as in real space. The presence 
of a surface breaks translational invariance, so that both 
the order parameter and the density distribution become 
dependent on position close to the surface. In partic- 
ular they exhibit "Friedel-like" oscillations. In Figs. 1 
and 2, this is demonstrated for the order parameter for 
a 64-site system with OBC, where site 1 and 64 repre- 
sent the surfaces. As mentioned above, the only distinc- 
tion in our tight-binding formalism is that there is no 
hopping matrix element connecting the surface sites to 
outside of the chain. Further refinements, such as alter- 
ing the hopping matrix elements or the attractive inter- 
action for sites near the surface, are possible; however 
we do not explore these possibilities here. Unlike the 
Ginzburg-Landau approaches, we do not impose bound- 
ary conditions on the order parameter itself: the surface 
behaviour of the order parameter is determined by the 
BdG equations themselves. 

In Figs. 1 and 2, the BdG and Anderson results are 
shown with solid and dashed curves, respectively. Fig- 
ure 1 shows the order parameter Aj (scaled by the site- 
averaged value) as a function of site number i for half 
filling, for various coupling strengths. At half filling, the 
order parameter is peaked at the surface and oscillates 
as one moves in towards the bulk. These oscillations de- 
cay over the scale roughly determined by the coherence 
length £, beyond which the order parameter relaxes to its 
bulk value. For sufficiently large samples, the bulk value 
is the same as that achieved for the homogeneous case 
obtained with PBC. This can be seen clearly in Fig. 1, 
with £ becoming longer as the coupling strength is de- 
creased. An exception is the case with \U\ = Lit, where 
we encounter finite size effects - the system is not quite 
large enough to allow relaxation of the order parameter 
to its bulk value. 

The "Friedel-like" oscillations in the order parameter 
are a reflection of the single-particle wave function at the 
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Fermi level. For a given density n, the oscillation pe- 
riod in site number is given by ir/kpa, where a is the 
lattice spacing, and &f and n are connected by the usual 
Pauli consideration. To illustrate this, we plot the order 
parameter with \U\ = l.lt for various values of density 
n in Fig. 2. As n decreases, Uf decreases and the pe- 
riod becomes longer, as is indicated by the bar in each 
panel. Note moreover that the coherence length £, for 
the given coupling strength \U\ = l.lt, decreases mono- 
tonically from n = 1.0 to 0.25: this is shown in Fig. 3(a), 
where we plot £ as a function of electron density n. In 
Fig. 3(b) £ for half filling is plotted as a function of cou- 
pling strength \ U\. In these figures, we have also included 
some analytical forms valid in the weak-coupling limit 
(£bcs), and in (b), one in the strong-coupling limit as 
well. Note, particularly in (b) the excellent agreement 
of the analytical forms with £ in the two extremes. It 
is clear that in both Figures 1 and 2 the Anderson pre- 
scription (dashed curves) captures the essential features 
of the "Friedel-like" oscillations in the order parameter. 
This is also the case for the electron density distribution 
(not shown). 

In Fig. 4 we examine the local density of states calcu- 
lated with both the BdG and Anderson formalisms. The 
delta functions in Eq. (7) have each been replaced by a 
normalized Gaussian, and the smoothing width for the 
results shown in all the LDOS figures is O.li. It can be 
seen that the Anderson prescription reproduces the BdG 
results remarkably well in all details. In this figure the 
LDOS at various sites is shown for the same system as in 
Fig. 1 but with \U\/t = 2.0. Near a surface the LDOS is 
quite different from that in the bulk (as represented by 
the LDOS in the middle of the sample, or by the aver- 
age DOS over all sites). Although the overall LDOS as a 
function of energy is very different from the bulk DOS, 
the energy gap in the spectrum on the surface is almost 
the same as that in the bulk (compare the BdG results 
for site 1 and 32). On the other hand, it is intriguing that 
the energy gap is larger at every second site (site 2 and 
4 in Fig. 4): this statement applies only to half filling. 
For quarter filling, for example, the energy gap is slightly 
larger at every fourth site. The Anderson results repro- 
duce these features correctly, and their agreement with 
the BdG results becomes better further away from the 
surface. Note, however, the absence of coherence peaks 
in the more accurate BdG results at the surface position. 

B. Friedel-like oscillations 

Next we study the "Friedel-like" oscillations in more 
detail. In Fig. 5 the density distribution func- 
tion of site number i is shown for a 64-site chain (OBC, 
\U\/t = 1.5) for various values of the average electron 
density n. In this figure the density distribution obtained 
from the BdG results (solid curves) is compared with the 



density distribution for the noninteracting case (dotted 
curves). Except for half filling, the density distributions 
in both the superconducting and normal states exhibit 
oscillations, with period determined by kp 1 . In the su- 
perconducting state the oscillations are smoothed as one 
moves away from the surface, as was the case with the or- 
der parameter (see Fig. 2). This distance over which the 
electron density becomes smooth is the same as for the 
order parameter, and is loosely given by the coherence 
length £. The exceptional case of half filling (n = 1.0) 
has an electron density that is independent of site posi- 
tion, due to the particle-hole symmetry present at half 
filling (see next subsection, however). 

As seen in Fig. 5, the density distribution (for nonzero 
U) sometimes has higher amplitude near a surface, par- 
ticularly close to half filling. This behaviour is not to 
be confused with the decay of the "Friedel-like" oscilla- 
tions in the order parameter. The former simply follows 
the non-interacting distribution, which results from the 
particle-in-a-box problem subject to the Pauli exclusion 
principle. With OBC the single-particle states are stand- 
ing waves, and the sum of their probability distribution 
tends to be higher near a surface near half filling. 

While the electron-electron interaction is local in 
this model, in two and three dimensional systems, the 
"Friedel-like" oscillations appear in somewhat different 
ways from one dimension (ID). The difference stems 
mainly from two factors. One is that unlike in ID, with 
OBC, most of the single- particle energy levels are degen- 
erate in higher dimensions. As a result, interference oc- 
curs among degenerate single-particle states at the Fermi 
level, and the order parameter and density distribution 
can exhibit quite complicated structure. In two dimen- 
sions (2D) at half filling, however, it is relatively simple 
due to the structure of the Fermi surface. We illustrate 
this in Fig. 6, where the order parameter is plotted for 
a 20 x 20-sitc system for \U\/t = 2.5 and 1.5. At half 
filling, the order parameter shows the "Friedel-like" os- 
cillations, where it is peaked along the diagonal sites and 
largest near the corners. As a function of decreasing \U\ 
(increasing £), the oscillations extend towards the middle 
of the sample. These oscillations, as seen in Fig. 6, can 
be explained by the interference of single-particle wave 
functions. In the zero-coupling limit (\U\ — > 0), Eq. (5) 
reduces to A 4 ~ \U\ (1/2) J]k F u k F («)v kF M at half filling 
at zero temperature, and uu F (i) — Vk F {i) — u k F M- Here 
u^ p (i) is the single-particle wave function (analytical so- 
lution [5]) for the Fermi momentum kp — (k x ,k y ), and 
the ~ denotes that the oscillatory behaviour is similar. 
In fact the amplitude of the hole-like part decreases with 
\U\ for electron-like states and vice- versa for the electron- 
like part.. Using the relation k y = — k x + tt on the Fermi 
surface to eliminate k y , and approximating the sum by an 
integral over k x , we can calculate Aj analytically. Thus 
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A,; ~ A 



1 + ik<W< + s vi,- 



(8) 



where x^, are the x and y coordinates of site i, and A is 
a small number. Clearly Aj is larger along the diagonals 
(xi = ±j/»). We have obtained the BdG solutions for 
N = 32 2 at half filling as a function of \U\, and found 
that for very weak coupling (e.g., \U\/t — 0.1), the above 
equation indeed describes the order parameter structure. 

For three dimensional systems, the largest size we have 
examined is TV = 12 3 . For coupling strengths small 
enough to see the "Friedel-like" oscillations, however, this 
size is still small - there are only 12 sites on one side and 
the surfaces occupy a large portion of the system - so 
that the order parameter is dominated by finite size ef- 
fects. We have calculated the three-dimensional version 
of Eq. (8) numerically, and found that it qualitatively ex- 
plains the oscillatory structure of the order parameter in 
the weak-coupling limit. 

As n becomes smaller than 1, the interference pattern 
changes. As seen above, in 2D, for n ~ 1, the Fermi 
surface is approximately k y = —k x + n and the oscilla- 
tions in Aj as well as n, are along several sites around 
the diagonal. Away from half filling, the oscillations oc- 
cur everywhere along the surfaces. Figure 7 shows the 
order parameter and density distribution for N = 20 2 
and \U\/t — 2 at quarter filling. In the low-density limit, 
both rii and Aj are suppressed on and near the surfaces. 

Another factor that distinguishes one dimension from 
higher dimensions is the dependency of nesting on dimen- 
sionality. In ID, the Fermi surface is perfectly nested for 
all values of n. In 2D and 3D this is the case only at 
half filling, and this means that away from half filling, 
"Friedel-like" oscillations can be suppressed compared 
with ID. We have not identified any noticeable differ- 
ence in this sense between ID and 2D. For 3D, we have 
examined this aspect by placing a nonmagnetic impurity 
which shifts the chemical potential at that site and using 
PBC. In this way a disturbance that breaks translational 
invariance is one site, as opposed to all the surfaces. In 
Fig. 8 we show the order parameter for N — ll 3 , n = 0.9, 
and \ U\/t = 2, in the first layer (top) and the second layer 
(bottom). The impurity potential is 0.5t in the middle of 
the first layer. The basic features caused by an impurity 
are the same as in ID [5]: in the first layer Aj is sup- 
pressed at the impurity site and shows the "Friedel-like" 
oscillations around it. The size-dependent electron den- 
sity, rij (not shown), behaves similarly, except that for 
an attractive potential, it is peaked at the impurity site. 
The Aj is peaked right below the impurity site in the sec- 
ond layer and suppressed in the third layer, and so on. 
However, the extent to which the oscillations propagate 
is smaller compared with ID and 2D. This is the case 
within the layer where the impurity is, as seen in Fig. 8. 
Moreover, the oscillation amplitudes decrease quickly as 
one moves away from that layer: note that in Fig. 8 the 



z-axis range for the second layer is half of that for the 
first layer. 

In the strong-coupling limit, at half filling in any di- 
mension, the order parameter is larger on the surfaces 
(largest at the corners) and more or less flat inside. In 
Fig. 9 we show A 4 for N = 12 3 and \U\/t = 2 at half 
filling, for the first four layers. While Aj is larger on the 
surfaces (the first layer and peripherals of each layer), 
some oscillations near the surfaces can be seen. 



C. CDW 

For the half-filled attractive Hubbard model with near- 
est neighbour hopping on a translationally invariant sys- 
tem, the superconducting (SC) and charge-density-wave 
(CDW) states are degenerate. The presence of an im- 
purity tends to make a CDW state more favourable, and 
one might expect that a surface has the same effect. This 
is not necessarily the case, however. With an even num- 
ber of sites N, with OBC, the SC and CDW states are 
degenerate at half filling. Moreover, the state with con- 
stant density distribution discussed above is not the only 
SC solution. When the density distribution of a CDW 
state is used as input, the BdG equations converge to 
a SC state with a similar density distribution, and this 
state has the same energy as the one with constant den- 
sity. The order parameter structure is the same for both 
states, but the overall magnitude is smaller for the for- 
mer. For any average density, the BdG equations always 
converge to a SC solution. 

The situation is different, however, for an odd-iV sys- 
tem with surfaces. The BdG equations converge to a SC 
solution at any density that corresponds to an odd num- 
ber of electrons, e.g., at half filling. For an even number 
of electrons, the BdG equations either do not converge at 
all, or converge to a CDW solution. This presumably is 
caused by the incommensurability of the electron num- 
ber (all paired) and the number of sites. The same CDW 
solution for an even number of electrons is obtained by 
solving the Hartree equations (no superconducting order 
parameter) . With N + 1 electrons the system is equiva- 
lent to an N + 1 system with PBC at half filling, with a 
strong repulsive impurity at one site where the electron 
density becomes zero. While the Hartree equations usu- 
ally do not converge for an odd electron number, we were 
able to find a CDW solution at half filling with odd N. It 
is interesting to note that unlike even N, the CDW state 
has a lower energy than the SC state, and we did not find 
a SC state with an oscillating density distribution. 

In Fig. 10(a) we show the average density n vs. chem- 
ical potential ("dressed" with the Hartree energy), ob- 
tained by BdG for ID with N = 15 with \U\/t = 1. The 
plateaus correspond to CDW states, where the supercon- 
ducting order parameter is zero. For stronger coupling 
and larger system size, the plateaus tend to disappear, 
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but the BdG equations pick up the CDW solutions at cer- 
tain densities. In Fig. 10(b) the density distribution is 
shown for a CDW state in ID with N = 63, \U\/t = 0.5, 
and n = 18/63. In both (a) and (b), the BdG and An- 
derson results are plotted with solid and dashed curves, 
respectively. Although the Anderson results are similar 
to the BdG results, they are not CDW solutions: the 
density distribution is the same as the one in the nonin- 
teracting case. 

D. Local density of states 

To distinguish the features that are due to the open- 
ing of a superconducting gap from those that are due 
to the single-particle states, we compare, in Fig. 11, the 
LDOS (the same BdG results as shown in Fig. 4) with 
that in the noninteracting case, for various sites near the 
surface. Interestingly, the noninteracting LDOS has a 
gap- like feature at every second site (e.g., site 2 and 4) 
near a surface (again this is the case only for half filling) . 
This is consistent with a larger gap at these sites in the 
superconducting state. 

As seen in Figs. 4 and 11, the energy gap as it appears 
in the LDOS on the surface is not very different from 
the one in the bulk, despite the fact that both the or- 
der parameter and the density distribution behave quite 
differently on and near the surface. However, the BCS 
coherence factors, which give rise to peaks in the density 
of states just beyond the gap (as seen in the LDOS for 
site 32 in Fig. 11), are barely existent near the surface; 
for half filling, particularly at the even-numbered sites. 
For an orbital structure much more complicated than the 
single-band one-dimensional model adopted here, these 
effects would have to be considered when interpreting 
STM or ARPES data. 

The site variation of the LDOS also depends on elec- 
tron density. In Fig. 12 we plot the LDOS at the surface 
of a 64-site chain for various average electron densities. 
We have solved the BdG equations self-consistently for 
\U\/t = 2.0 and for the given values of electron density. 
The resulting LDOS is compared to that calculated in 
the bulk (we simply averaged over all sites, which gives 
a result identical to that obtained in the interior, or that 
obtained with PBC) . Particularly in the low-density limit 
large differences occur: the bulk DOS always has a well- 
defined energy gap, with remnants of the BCS coher- 
ence peaks, whereas the LDOS on the surface does not. 
This occurs mainly because in the low-density limit, the 
electron density at the surface sites is significantly lower 
than in the bulk. This is most apparent in Fig. 5, where 
the electron density (even in the non-interacting limit) 
is considerably lower near the surface than in the bulk. 
This is true even when we include repulsive Coulomb 
interactions between electrons, as we have verified by di- 
rect calculation in the low-density limit. In other words, 



Coulomb interactions do not affect the inhomogeneous 
electron distribution as much as one might think because 
we are in the dilute limit. 

In 2D and 3D, the LDOS behaves in a similar way as in 
ID, as one moves away from a corner along the diagonal. 
In Fig. 13 the LDOS for N = 32 2 and \U\/t = 2.5 at 
quarter filling is shown for several sites near the surfaces. 

IV. SUMMARY 

We have calculated the spatially inhomogeneous or- 
der parameter and local density of states (LDOS) for a 
very simple model exhibiting s-wave superconductivity, 
as a function of coupling strength and electron density. 
We have used the Bogoliubov-de Gennes formalism and 
solved the resulting equations self-consistently, including 
both the anomalous and normal effective potentials. We 
find that both the order parameter and the density dis- 
tribution vary significantly as a function of position near 
the surface. The details of this variation are dependent 
on both the electron density and the coupling strength. 
The LDOS shows a similar variation. If one can some- 
how preferentially tunnel electrons into the surface or 
subsurface layers, this variation may be observable [14]. 
Here again, details will be dependent on the specifics of 
the model, i.e., the character of the orbitals involved in 
the superconductivity as well as the symmetry of the su- 
perconducting order parameter. We intend to further 
investigate these issues in future work. 
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FIG. 1. Order parameter A, normalised by the average 
value as a function of site number i, for a 64-site chain with 
OBC at half filling, for various values of coupling strength 
\U\. The BdG and Anderson results are plotted with solid and 
dashed curves, respectively, which agree very well. The gap 
exhibits the "Friedel-like" oscillations near a surface, which 
decay roughly over the coherence length £ (indicated by ar- 
rows). 
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FIG. 2. Same as Fig. 1, but for \U\ = Lit and for various 
electron densities n. The period of the "Friedel-like" oscilla- 
tions in site number i is roughly n/kpa, where a is the lattice 
constant, and increases as n is decreased. The irregular oscil- 
lations seen for n = 0.75 arise from the incommensurability 
of the electron number and the number of lattice sites. The 
coherence length £ (see next figure) decreases as a function of 
decreasing n. 
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FIG. 3. The BCS coherence length as a function of (a) 
electron density, and (b) coupling strength, with comparison 
to analytical approximations. 
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FIG. 4. The LDOS for an N = 64 system with OBC at 
half filling, with \U\/t = 2.0, for several sites: site 1 is a 
surface and site 32 is the middle of the sample. The Anderson 
approach (dashed curves) reproduces the BdG results (solid 
curves) remarkably well. Near a surface, the LDOS is quite 
different from that in the bulk (site 32), and at half filling, 
the energy gap is larger at every second site (e.g., site 2 and 
4). 



FIG. 5. Density distribution m for an N = 64 system with 
OBC and \U\/t = 1.5 for various average electron densities 
n. The BdG results (solid curves) are compared with those 
for the noninteracting case (dotted curves). Except for half 
filling, rii shows "Friedel-like" oscillations near a surface, fol- 
lowing the noninteracting density distribution. Similarly to 
A; , the length scale for the decay of these oscillations is given 
by the coherence length £ (indicated by arrows). 




FIG. 6. Order parameter Aj, as a function of position on a 
2D (20X20) lattice with OBC, for \U\/t = 2.5 (upper graph) 
and \U\/t = 1.5 (lower graph). Note the variation near the 
surfaces and along the diagonals, particularly for the weaker 
coupling case. These results are computed at half filling. 
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FIG. 7. Order parameter A; (upper graph) and n(i) (lower 
graph), as a function of position on a 2D (20 x 20) lattice with 
OBC, for \U\/t = 2. These results are computed at quarter 
filling (n = 0.5). Note the variation near the surfaces. 



FIG. 9. Order parameter Aj, as a function of position on 
a 3D (12 x 12 x 12) lattice with OBC, for \U\/t = 2, at half 
filling, n = 1. Results are shown for four layers, starting with 
the surface layer. The order parameter displays considerable 
variation near all the surfaces. 
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FIG. 8. Order parameter A», as a function of position on 
a 3D (11 x 11 x 11) lattice with PBC, for \U\/t = 2, with the 
average electron density, n = 0.9. A single-site impurity is 
located in the middle of the surface layer (upper graph) with 
strength e = 0.5t. The variation of the order parameter is also 
shown in the second layer (lower graph) . Note the condensed 
scale in the lower graph, indicating that the disruption of the 
order parameter is concentrated near the impurity. 



FIG. 10. (a) Average electron density vs. chemical poten- 
tial and (b) electron density as a function of position, for 
an odd number of sites N with OBC. Solid curves (dashed 
curves) are BdG (Anderson) results. We used \U\ = It and 
N = 15 in (a) and \U\ = 0.5*, N = 63, and n = 18/63 in (b). 
The latter graph is for a CDW state (with zero superconduct- 
ing order parameter). 
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FIG. 11. Same as in Fig. 4, but the BdG results (solid FIG. 13. Site dependent LDOS for N = 32 2 and 

curves) are compared with those for the noninteracting case \U\/t = 2.5 at quarter filling (n = 0.5), for several sites near 
(dotted curves). the corner (along the diagonal). The dot-dashed curves are 

the average (bulk) result. 




FIG. 12. The LDOS for an N = 64 chain with OBC and 
\U\/t = 2.0, for several average electron densities n, obtained 
by self-consistent solution to the BdG equations. Those at 
the surface (solid curves) are compared with the average for 
all sites (dash-dotted curves), which has the same behaviour 
as the bulk DOS. The former differ markedly from the latter, 
particularly in the low-density limit. 
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